/* julia set fractals */

#include <stdio.h>
#include <stdlib.h>
#include "print_png.h"
//#include <cilk-lib.cilkh>
#include <string.h>
#include <complex.h>
#include <math.h>
#include <pthread.h>

// Maximum number of iterations
#define MAX_ITER 255
#define THRESHOLD 10

#define POLY 1
#define SINE 2
#define COSINE 3

#define MAXTHREADS 32

typedef struct t_args{
  double xMin;
  double xMax;
  double yMin;
  double yMax;
  int size;
  complex double complexIn;
  int type;
  int start;
  int finish;
  double *coef;
  char *values;
}t_args_t;

extern double wtime();

/* Function calc_mb
   Calculates the number of iterations to determine whether or not a point
   is a part of the mandlebrot set
   Args xMin, xMax, yMin, yMax - map out the area for which the mandlebrot
     set is being calculated.  Not changed.
   Arg values - array in which the calculated values will be stored.
     Values are overwritten on output.
   Arg size - size of the nxn values array.  Not changed.
 */

void *t_calc_julia(void *ptr)
{
  printf("Creating thread\n");
  struct t_args *args = (struct t_args *)ptr;
  calc_julia(args->xMin, args->xMax, args->yMin, args->yMax,
	     args->values, args->size, args->complexIn,
	     args->type, args->coef, args->start,
	     args->finish);
}


int calc_julia(double xMin, double xMax, double yMin,
		    double yMax, char *values, int size,
		    complex double complexIn, int type, double* coef,
		    int start, int finish){

  int i, j, k, result;
  double x, y;
  complex double point;
  complex double complexNum;

  printf("start: %d  finish : %d\n",start,finish);
  
  // if(finish-start<50){

    for(i=start; i<finish; i++){
      for(j=0; j<size; j++){
	x = ((double)i/size)*(xMax-xMin) + xMin;
	y = ((double)j/size)*(yMax-yMin) + yMin;
	point = x + y * _Complex_I;
	//point = Complex(x,y);
	//	point.r = x;
	//point.i = y;
	complexNum = complexIn;
	for(k=0; k<MAX_ITER; k++){

	  if(type==POLY){
	    point = coef[3] * cpow(point, 3.0) + coef[2] * cpow(point, 2.0) +
	      coef[1] * cpow(point, 1.0) + coef[0] + complexNum;
	    //	  complexNum = coef[3] * cpow(complexNum, 3.0) +
	    //  coef[2] * cpow(complexNum, 2.0) + coef[1] * complexNum + 
	    //  coef[0] + point;
	  }
	  if(type==SINE){
	    point=complexNum*csin(point);
	    // complexNum = point*csin(complexNum);
	  }
	  if(type==COSINE){
	    point = complexNum*ccos(point);
	    //complexNum = point*ccos(complexNum);
	  }
	  // rule of thumb cutoff that shows a point will not be in the set
	  if(cabs(point)>THRESHOLD){
	    //	  printf("cabs = %f\n", cabs(complexNm));
	    break;
	  }
	  
	}
	values[j*size+i]=k;
      }
    }
    return 0;
}


int main (int argc, char* argv[]){
  
  // Will be used to create the image
  char *values; 

  // Number of pixels in the generated square
  int size;
 
  // default values
  double xMin = -1.5;
  double xMax = 0.5;
  double yMin = -1.0;
  double yMax = 1.0;

  double t1, t2, diff;

  char fileNameBase[50];
  char fileName[50];
  char line[20];

  double coef[4];
  int i, result, start, finish, blockSize;

  // input file
  FILE * infile;

  // complex
  double real;// = -0.123;
  double imag;// = 0.745;
  complex double complexNum;

  struct t_args thread_args[32];

  pthread_t tid[MAXTHREADS];

  if((argc!=4)){
    printf("Usage: julia.x $SIZE $OUTFILENAME $INFILE\n");
    exit(1);
  }

  infile = fopen(argv[3], "r");
  if(infile==NULL){
    printf("Error opening input file\n");
    exit(1);
  }
  else{
    printf("Reading from %s\n",argv[3]);
  }

  // complex point to calculate set for
  // fuck knows why i can't get fscanf to work right
  //fscanf(infile, "%f", &real);
  fgets(line, 20, infile);
  real = atof(line);
  fgets(line, 20, infile);
  //  fscanf(infile, "%f", &imag);
  imag = atof(line);
  complexNum = real + imag * _Complex_I;
  //complexNum = Complex(real,imag);
  //complexNum.r = real;
  //complexNum.i = imag;

  // box
  fgets(line, 20, infile);
  xMin = atof(line);
  fgets(line, 20, infile);
  xMax = atof(line);
  fgets(line, 20, infile);
  yMin = atof(line);
  fgets(line, 20, infile);
  yMax = atof(line);

  // polynomial
  for(i=3;i>=0;i--){
    fgets(line, 20, infile);
    //    fscanf(infile, "%f", &temp);
    coef[i]=atof(line);
  }

  fclose(infile);
  
  if(1){
    printf("Read in: \n");("Point = %f + %fi\n", creal(complexNum),cimag(complexNum));
    printf("Box: xvals %f-%f, yvals %f-%f\n",xMin,xMax,yMin,yMax);
    printf("Poly: %fx^3 + %fx^2 + %fx + %f\n",coef[3],coef[2],coef[1],coef[0]);
  }

  size=atoi(argv[1]);
  
  strcpy(fileNameBase,argv[2]);

  values = (char*) malloc (sizeof(char)*size*size);
  if(values==NULL){
    fprintf(stderr,"Malloc failed\n");
    exit(1);
  }

  t1=wtime();

  blockSize=size/(MAXTHREADS-1);

  for(i = 0; i<MAXTHREADS; i++){
    thread_args[i].xMin=xMin;
    thread_args[i].xMax=xMax;
    thread_args[i].yMin=yMin;
    thread_args[i].yMax=yMax;
    thread_args[i].values=values;
    thread_args[i].size=size;
    thread_args[i].complexIn=complexNum;
    thread_args[i].coef=coef;
    
    thread_args[i].type=POLY;
  }
  
  for(i=1;i<MAXTHREADS;i++){
    start=(i-1)*blockSize;
    finish=i*blockSize;
    
    thread_args[i].start=start;
    thread_args[i].finish=finish;
    
    pthread_create(&tid[i], NULL, (void *(*)(void *))t_calc_julia, (void *)&thread_args[i]);
  }
  //calc_julia(xMin, xMax, yMin, yMax, values, size, complexNum, POLY, coef, 0, size);

  for(i=1;i<MAXTHREADS;i++){
    pthread_join(tid[i],NULL);
  }

  strcpy(fileName, "POLY");
  strcat(fileName, fileNameBase);
  print_png(values, size, size, fileName);


  for(i=1;i<MAXTHREADS;i++){
    start=(i-1)*blockSize;
    finish=i*blockSize;
    thread_args[i].type=SINE;
    thread_args[i].start=start;
    thread_args[i].finish=finish;
    pthread_create(&tid[i], NULL, (void *(*)(void *))t_calc_julia, (void *)&thread_args[i]);
  }

  for(i=1;i<MAXTHREADS;i++){
    pthread_join(tid[i],NULL);
  }

  //  calc_julia(xMin, xMax, yMin, yMax, values, size, complexNum, SINE, coef, 0, size);
  
  strcpy(fileName, "SINE");
  strcat(fileName, fileNameBase);
  print_png(values, size, size, fileName);

  for(i=1;i<MAXTHREADS;i++){
    start=(i-1)*blockSize;
    finish=i*blockSize;
    thread_args[i].type=COSINE;
    thread_args[i].start=start;
    thread_args[i].finish=finish;
    pthread_create(&tid[i], NULL, (void *(*)(void *))t_calc_julia, (void *)&thread_args[i]);
  }

  for(i=1;i<MAXTHREADS;i++){
    pthread_join(tid[i],NULL);
  }

  //  calc_julia(xMin, xMax, yMin, yMax, values, size, complexNum, COSINE, coef, 0, size);
  strcpy(fileName, "COSINE");
  strcat(fileName, fileNameBase);
  print_png(values, size, size, fileName);
  
  t2=wtime();



  diff=t2-t1;

  printf("Time to run program = %f\n",diff);
  printf("%dx%d matrix\n",size,size);
  printf("xmin = %f, xmax = %f\n",xMin,xMax);
  printf("ymin = %f, ymax = %f\n",yMin,yMax);

  return 0;
}

